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Abstract 



We study associative memory based on temporal coding in which successful 
retrieval is realized as an entrainment in a network of simple phase oscilla- 
tors with distributed natural frequencies under the influence of white noise. 
The memory patterns are assumed to be given by uniformly distributed ran- 
dom numbers on [0, 2tt) so that the patterns encode the phase differences of 
the oscillators. To derive the macroscopic order parameter equations for the 
network with an extensive number of stored patterns, we introduce the ef- 
fective transfer function by assuming the fixed-point equation of the form of 
the TAP equation, which describes the time-averaged output as a function of 
the effective time-averaged local field. Properties of the networks associated 
with synchronization phenomena for a discrete symmetric natural frequency 
distribution with three frequency components are studied based on the order 
parameter equations, and are shown to be in good agreement with the results 
of numerical simulations. Two types of retrieval states are found to occur 
with respect to the degree of synchronization, when the size of the width of 
the natural frequency distribution is changed. 
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I. INTRODUCTION 



Recently a population of neurons in the cat visual cortex has been reported to exhibit 
synchronized firings in a stimulus dependent manner . The occurrence of correlations in 
firing times of neurons seems to be a ubiquitous phenomenon in real nervous systems. The 
role of such synchronized firings for information processing in the brain has been attracting 
growing interest of researchers, and several authors have suggested neural network models 
based on the concept of temporal coding, where information of a neuron is represented by 
its firing times. Indeed, to explain the experimental findings that visual information of an 
external object is processed with being divided into several pieces of information, several 
authors [|]-|(| have suggested that the synchronized firings of neurons may serve as a linker 
for those pieces of information. 

The problem of investigating how an associative memory is realized in real nervous sys- 
tems as well as of constructing biologically relevant models is of central concern of neuroscien- 
tists. Since the establishment of systematic theories of associative memory for networks with 
an energy function that is ensured by assuming symmetric synaptic couplings f?HT6l , several 
attempts have been made to make models as biologically plausible as possible [p| -|T9|j . Pre- 
viously we investigated the effects of asymmetric couplings [^Of for memorizing pre-synaptic 
and post-synaptic activities, which are incorporated into the standard symmetric Hebb 
learning rule [21|, by studying networks of analog neurons p2,[23|, whose continuous-time 



dynamics involves a positive- valued transfer function representing the mean firing rates of a 
neuron as a function of membrane potential. 

Instead of working with the concept of rate coding based on the idea that neuronal 
information is represented by mean firing rate of a single neuron, one may be concerned with 
the concept of temporal coding, when considering that spatio-temporal patterns of neuronal 
firings will make information carried by a population of neurons much richer than spatial 
patterns alone. A spiking neuron is considered to be one of the candidates for implementation 
of temporal coding p4|-|33|j. The time evolution of membrane potential that is generated in 
response to an injected synaptic electric current of a spiking neuron is described by such 



nonlinear dynamics as Hodgkin-Huxley equation FitzHugh-Nagumo equation P5] , |56| , 
or the equation of an integrated-and-fire neuron. Spiking neurons in a network are generally 
supposed to interact with each other via pulses generated in the firing events occurring in 
the pre-synaptic neuron. We have shown previously that even in the presence of time delays 
in transmission of the pulses an associative memory based on a network of spiking neurons 
can be realized by assuming a simple Hebb-type learning rule alone, and that the memory 
retrieval accompanies synchronized firings of neurons. The dynamics of such associative 
memory was analyzed by means of sublattice method in our previous paper [37] . 

In the previous analysis |37| we assumed that every neuron shares identical characteristics 
to exhibit the same reaction in response to the same injected current. In real nervous 
systems, however, neurons in a network may have their own individual characteristics. 

The problem of whether temporal coding functions robustly in the presence of certain 
heterogeneties of neurons will be of particular interest. For the purpose of investigating such 
a problem for the phenomenon of synchronized oscillations in associative memory neural 
networks, we consider it appropriate to deal with simple models of coupled oscillators with 
distributed natural frequencies and external noise. It is well known that, under certain 
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conditions, a population of oscillators with distributed natural frequencies is allowed to get 
partially entrained in such a way that oscillators with a natural frequency near the central 
frequency become to oscillate synchronously at identical frequency as a result of cooperative 
interactions [p^ -|40| • 

Kuramoto |38| showed that the dynamics of this kind of network of oscillators with 
sufficiently week interactions can be reduced to a simple phase dynamics. Supposing that 
neurons in a network are treated as phase oscillators, associative memory has been shown 
to be realized under a simple learning rule of the Hebb-type either in the case of a finite 
number of stored patterns [£I,39| or in the case with a single natural frequency p2,£5 



Satisfactory analysis of the case with a distribution of natural frequency and extensively 
many stored patterns has been far less conducted. Quite recently we have reported the 
study of a deterministic phase oscillator network with a distribution of natural frequencies 
where an extensive number of binary patterns (±1) are stored with use of the Hebb learning 
rule 0. 

The main purpose of the present study is the theoretical analysis of associative memory 
based on temporal coding with use of networks of phase oscillators in the more general case 
where the number of stored patterns that are given by uniformly distributed random numbers 
on [0, 27r) is extensive, natural frequencies of the oscillators are distributed according to a 
certain distribution function, and furthermore external white noise is added to the system. 

While one can analyze a phase oscillator network with a single natural frequency by 
means of the replica method that makes full use of its associated energy function, one 
cannot resort any more to the standard method of statistical physics based on the existence 
of an energy function in the case of networks with a distribution of natural frequencies. 

One can, however, use the Self-Consistent Signal-to-Noise Analysis (SCSNA) 
^5| , |^B|P^ , |25|J3T| to deal with general cases without energy functions. To apply the SC- 



SNA it is necessary to know fixed-point equations describing the equilibrium states of the 
network. When considering such equations in stochastic systems, we may take advantage of 
the concept of the TAP equation [£7|, 31l| . 



The TAP equation is known to exist for the SK model of spin glasses j48H50j and the 
Hopfield model of an Ising spin neural networks B, f48| , |5T|| , and to represent a functional rela- 
tionship between the thermal or time average of each spin in equilibrium and its correspond- 
ing effective local field that involves the so called Onsager reaction term [f48[]. Usefulness of 
the TAP equation in deriving the order parameter equations of associative memory networks 
is attributed to the fact that the resulting equations of the replica calculations by AGS is 
recovered by the result of application of the SCSNA to the TAP equation where the 
Onsager reaction term is canceled exactly by the appearance of the renormalized output 
term of the SCSNA | 46|| . We note that the TAP equation of the naive mean field model 
with the interactions given by the Hebb learning rule defines to an analog network equation 
with the transfer function tanh (j3h) [ [46|j . 

We first evaluate an analogue of the TAP equation of the naive mean field type for our 
model with a distribution of natural frequencies by dealing with the Fokker-Plank equation. 
In order to obtain the Onsager reaction term we compute the free energy of the network 
without a distribution of natural frequencies to derive the TAP equation by following the 
method of Plefka |50| and Nakanishi |5T| . Then we assume that the TAP-like equation 



also exists with the Onsager reaction term remaining the same even for networks with 
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a distribution of natural frequencies, and that such TAP-like equation defines an effective 
transfer function to which the SCSNA is applicable to obtain the order parameter equations. 

The present paper is organized as follows. In section |I| we introduce a neural network of 
simple phase oscillators and describe how the network functions as an associative memory 
based on a simple learning rule of the Hebb-type. In section [TT1] we give a theoretical analysis 
to derive the macroscopic order parameter equations describing the long time behavior of 
the system. On the basis of the order parameter equations, in section [TV] we investigate 
properties of memory retrieval accompanying synchronization in the networks by assum- 
ing a discrete symmetric natural frequency distribution with three frequency components. 
Results of numerical simulations are presented showing good agreement with those of theo- 
retical analysis. In section [V], comparing our work with those of other researchers conducted 
previously, we summarize the results of the present study. 

II. NEURAL NETWORKS OF PHASE OSCILLATORS WITH DISTRIBUTED 

NATURAL FREQUENCIES 

The system under consideration is a network of N phase oscillators subjected to external 
white noise, whose dynamics is expressed as 

N 

(pi = uJi-Y^ Jij sin (4>i - 4>j - Pij) + T]i(t), (1) 

where fa and Ut represent the phase and the natural frequency of oscillator i respectively. 
(3ij and represent a certain phase shift and the strength of coupling between oscillator 
i and j respectively. The Gaussian white noise rji(t) is assumed to satisfy {r)i(t)) = and 

When uji = (i = 1, . . . , N), = exp satisfy c y - = Cj* with * denoting complex 
conjugation and Jy take real values, the system (U) has the energy function: 

= -\ E Jn cos - h ~ /%)■ (2) 

Then one has an equilibrium probability distribution p{{fa}) proportional to 
exp (— H({<f>i})/D). In the case of D — 0, the function (@) becomes a Lyapunov func- 
tion of the system and hence the state of the system eventually settles into a certain fixed 
point attractor after a long time. 

In the present study we assume natural frequencies to be distributed accordingly to an 
even distribution function g{uj) = g(—uj) so that the average of natural frequencies become 
zero without loss of generality. To store P quenched random patterns 0f (i — 1, . . . , N, \i — 
1, . . . , P) chosen from the uniform distribution on the interval [0, 2tt), we assume the Hebb 
type learning rule, and set the parameter /3y and real valued such that 

Cij = J y exp (z/%) = I ^ «• ^ YJ., (3) 

where £f = exp (z^f). This definition of couplings gives networks the following properties: 
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• In successful retrieval an entrainment occurs, where synchronized oscillators satisfy 
(p i — <f)j pa 9} — 6j with a target pattern \x recalled. (Note that if fa is the solution 
of dynamics (0), uniformly shifted phase fa = (pi + c also becomes its solution. What 
matters is not the phase itself but their difference fa — <f>j.) 

• In the case of unsuccessful retrieval, all the oscillators fail to synchronize, running at 
their own natural frequencies. 

To measure the distance between the pattern \i and the state of the system, we introduce 
the overlap for pattern \i: 

= ^Etf**(*), ( 4 ) 

i 

where we denote by Zi(t). Then by use of the local field: 

Ht) = E^jW = Etf™^*) - ( 5 ) 

j+i n 

the dynamics (fj) is rewritten as 

fa = uji- Re{/ij(t)}sin0i + lm{hi(t)} cos fa + r)i(t), (6) 

where a denotes the loading rate P/N. From Eq. @ it is easy to see that the learning rule (||) 
indeed realizes the above mentioned properties if the number of stored patterns is finite (a = 
0) and Ui = (i = 1, . . . , N) . The memory retrieval accompanying synchronization can 
also occur for a > even in the presence of a distribution of natural frequencies. 

III. MACROSCOPIC ORDER PARAMETER EQUATIONS 

Behaviors of associative memory networks depend crucially on the nature of the local 
fields or the neurons of oscillators, because the updating rule for the time evolution of the 
system is based on the local fields as is seen in Eq. (^). When the long time behavior of 
a network is described by fixed point type attractors, the relation between the resulting 
output state of a neuron and the corresponding local fields becomes essential for determin- 
ing equilibrium properties of the networks. Such a relation is naturally introduced in the 
case of deterministic analog networks where neurons are characterized by transfer functions 
describing the input-output relation. 

For such stochastic systems ([I]) as Ising spin networks equilibrium fixed-point equations 
called the TAP equations are known to exist as expressing the relation between time average 
of each output of a neuron and its corresponding effective local fields involving the so called 
Onsager reaction term that is proportional to the time- averaged output. Once the TAP-like 
equation level description is available, the SCSNA, in which one compute the variance of 
the cross-talk noises in the local field as a result of storing an extensive number of patterns, 
can be applied to obtain the order parameter equations in the limit iV — ► oo. 
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A. Effective transfer function based on time-averaged local field and the TAP-like 

equation 



If Ui = (i = 1, . . . , N) and D = 0, the state of the network eventually settles into an 
equilibrium state given by a fixed point attractor owing to the existence of the Lyapunov 
function (Q), and then the local fields do not fluctuate in time. Even in the presence of 
external white noise (D ^ 0), the local fields also get fixed in time after a long time, 
provided a = P/N = 0. When the local fields are fixed over the time change due to the law 
of large numbers, theoretical treatment becomes simple because one can reduce the many 
body problem to a single body problem. 

In more general case where a > and/or natural frequencies are distributed, the local 
fields may fluctuate even with a large number of oscillators as can be shown in the numerical 
simulation illustrated in Fig. [I]. The fluctuations seem to be aperiodic and rigorous analysis 
of such fluctuations is quite difficult. To deal this situation we are forced to resort to a 
certain approximation by confining ourselves only to the near equilibrium behavior of the 
system: we replace the time- dependent local fields by their time-averaged ones: 

hi(t) ^f^^^Cijl-, (7) 

where the over bar represents the time average at near equilibrium. 

Once we apply this approximation, we are allowed to treat each interconnected oscillator 
as an element obeying the dynamics of one degree: 



<pi = uji- Re(h^ sin fa + lm(h^j cos fa + rji(t), (8) 

which is expected to describe the behavior of the oscillator i near equilibrium. It turns 
out that the time average Y[ at equilibrium can be expressed as a function of uji and hf 
~zi = f (uji, h{, h*^j . We call that the transfer function in the case of network of analog 

neurons. It is easy to show that the transfer function f(u, h, h*^j satisfies 

f{co,re i6 ,{re id Y) = e id f(u,r,r). (9) 
Hence it suffices to calculate f(u,r,r) to obtain f(ou,re ie , (re 10 ^ In the absence of white 



noise (D — 0), we can easily obtain the transfer function for real- valued local field r |38 



; f uj + v 'uj 2 — r 2 J uj < — r 
f(u,r,r)=^ ^ (ioj + -\A* 2 — w 2 ) —r<u<r . (10) 
: (u — y/uj 2 — r 2 ) r < uj 

In deriving Eq. ( |I0D in the case of —r < uo < r we used = together with d(j)/d(f) < 0. In 
the case of uo < —r or r < uj we computed the time average of exp (i<f>(t)) over the period 
T of the periodic oscillations of 0: 

/(", r, r) = ^ / exp (i0(t))dt = = (11) 

T ° Jo Jo jd<f> Jo 
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In Fig. ^|we illustrate the shape of the transfer function f(u, h, h*^j that is obtained from 
Eqs. (g) and 

In the presence of white noise (D > 0), from the Langevin equation (H) we obtain the 
Fokker-Plank equation: 

p((f),t; u,t) = -^{(u -Re(h)sin(P + lm(h)cos(p)p} + D^, (12) 

where p^0, t; oj, hj is the probability distribution of phase <p £ [0, 27r) at time t, and periodic 

boundary conditions p(o, t; u, hj = p(2rr, t; u>, JiJ, dp(o, t; uu, h\/d(f) = dp(2ir, t; u>, fij /d(j) are 
imposed. Since we are concerned with the probability distribution p eq (4>; u>, r) attained after 
a long time, we set p(0, t; ou, r) = to obtain |39| , |4l| 



p ^ r) = WJm^ (13) 

with 

I(4>) = exp (f cos 4>) / exp {— Hap — f cos Up + <p)}d^p, (14) 

JO 

where u = u/D and f = r / D. Noting the Ergodic property on the Fokker Plank Eq. (|T2"D, 
which holds when h is viewed as a given parameter, we obtain the time average of z by 
computing the average over the equilibrium distribution p eq (<f);uj,r): 

The transfer function we have obtained here can be considered to be an analogue of the so 
called TAP equation by the naive mean field model, because the time average of the output 
z is represented as a function of time-averaged local fields h. We note, however, that the 
genuine TAP equation, which is defined for systems with an energy function, should describe 
the functional relation between the time averaged output and the effective local field that 
differs in general from the time-averaged one. The difference between the two types of local 
fields is known to be the Onsager reaction term in the theory of random spin systems. We 
assume that the TAP-like equation may hold even for the present system without an energy 
function, and suppose it to be given by the following equations: 

z-=/(^,^ AP ,^ AP *), (16) 
^ AP =^ + 7 TAP ^ 

= C ii~3 + 7 TAP ^ 

= E^-«^ + 7 TAP ^. (17) 

A* 

In the case with uji = (i — 1, . . . , N) , by evaluating the free energy of the system, we can 
derive an explicit expression of the coefficient 7 TAP taking the form (see appendix |): 
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where 



?=^EN 2 (is) 

i 

It will, however, be difficult to rigorously derive an expression of 7 TAP for the general case 
with a natural frequency distribution. Thus we are led to make an assumption that the 
legitimate expression ( pT8[ ) for the case with u>i = (i = 1, . . . , N) can naturally be extended 
to the general case. Describing the desired 7 TAP requires the introduction of the order 
parameter u that appears in the SCSNA. 



B. Self-Consistent Signal-to-Noise Analysis 

We consider the case with m 1 = 0(1) and m M = 0(\/\fN^j (/i = 2, . . . ,P), where we 
choose pattern 1 as the target. Assuming, without loss of generality, that Q = 1 for all i 
and m 1 is real owing to the rotational symmetry, the local field Eq. ( |TTD is rewritten in the 
form: 

^ = ^ + ^EEerr^+7 TAP ^ (20) 

where we have used S{, hi, and m M to represent respectively zl, hi, and for brevity. 

When we consider the case with a finite number of stored patterns, the analysis is 
straightforward since we already know the form of the transfer functions ( |i0|) (D = 0) 
and fll~5|) (D / 0). Since 7 TAP = and hi = m 1 in this case, we have in the limit N — > oo 

ml = ^X#* a < - (/(wW)) , (21) 

where /. . .\ = J g{uj) . . .duj. Solving this equation numerically, the size of the overlap is 

evaluated as a function of various parameters including D. 

In the case of an extensive number of stored patterns (a > 0), however, cross-talk noise 
(the second term of Eq. (|2DD) in the local fields becomes to an appreciable extent, then 
one has to employ the method of SCSNA. The crux of the SCSNA is the evaluation of the 
variance of cross-talk noise that interferes occurrence of retrieval state. 



According to the prescription of the SCSNA p3| , the local fields ( P0|) is assumed to be 
in the form: 



hi = m 1 + E^r^ + 7 SCSNA ^ + 7 TAP *, (22) 



M>1 j¥=i 



where Sj, whose explicit expression is given later, is a quantity that is very close to and 
has negligible correlation in the limit N — > oo, and A and 7 SCSNA are to be self-consistently 
determined in the course of analysis. 



S 



Defining hi = m 1 + ^E M >i V; Eq. (pD reads 

* = f(u> it hi + 7 TOTAL ^, (hi + 7 TOTAL ^)> ( 23 ) 

where 7 TOTAL = 7 scsna + 7 tap_ 

Considering a general case with 7 TOTAL o, we solve Eq. (^) for s, to obtain the 
renormalized output Sj = f(uJi,h~i,hi J. Performing a Taylor expansion of f\0Ji,h~i,hi ) 

about (h~i,h~i ) , we have 



" £ JJ ^ (w) £ 



^Emr^l ( 24 ) 



with 



* = ^ + ^E Ear^- (25) 



N . . . 



ffaXXh (26) 



Substituting Eq. (|24]) into Eq. (p0|) and comparing the result with Eq. (|22|) (see appendix [IJ 
for details) we obtain 

^ + ^E E tftf**? + 7 SCSNA ^ + 7 TA % 
= m'+ E E & £f 4 + + 7 TA % (27) 



where 



1 ^5/ 



(28) 

hi, hi ) 



Since this equation holds for every site i, A and <-y b ° bJNA are self-consistently determined as 

A = rb (29) 

SCSNA = (30) 
1 — It 

The variance of renormalized noise iVj = ^X^>i Sjy* £,i£,j* s j m the right hand side of 
Eq. (p2|) can be evaluated by noting that RefiViJ and Im^/V^ distribute over sites obeying an 

identical Gaussian distribution independently and the site average of |iVj| 2 can be replaced 
by the pattern average: one has 

Tr T,\m 2 = (\Ni\ 2 ) =a\X\\ (31) 



N 



To represent the noise as N{ = ^Jarixi + iyi)/2, where Xi and yi are real and obey a normal 
Gaussian, we define 



r = 2|A| 2 g 



2q 



1 1 — u 



2 ■ 



(32) 



Summarizing Eqs. (|) , (|19|) , (II) > (HD ; © > (11) ; © , and fl32|), we have a set of macro- 
scopic order parameter equations 



s = f{co,h + 1 ^sAh + 1 



TOTAL 



h = m + 



x + iy), 



TOTAL SCSNA , „,TAP 

7 = 7 + 7 



m 



f(u,h,h* 



aru 



(x -iy)f(u,h,h* 
2q 



SCSNA 



where 



7 



represents 



1^? 



,2 » 



a- 



ff(^) 
2tt 



exp 



+ y 2 " 



(33) 

(34) 
(35) 
(36) 

(37) 
(38) 
(39) 
(40) 

. . . dudxdy. Detailed 



derivation of Eqs. ( p8|) and (|39|) is given in appendix pi . 

To discuss the generalized expression for 7 TAP for the case with a distribution of natural 
frequencies, we consider for the moment the case with ui = (i = 1, . . . , N) , where Eq. (]T8|) 
exactly holds. In this case it turns out that ^ TOTAL = by a rough argument given below. 
Note that ^ TOTAL = implies f (uj,h,h*^j = f(u,h,h*), and that the effective transfer 
function (JTJI) becomes 



/(0,r,r) 



Jo *" cos exp (f cos 
J 27r exp (f cos 



(41) 



Using Eqs. (|§P and flU]) and performing the average over the Gaussian distribution with unit 
variance for Eq. (|38|) , we obtain 



u= (l-q)/2D. 



Then, from Eq. we have 



7 



TAP 



U 



-oc- 



(42) 



(43) 



which immediately reconfirms 
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7 TOTAL = 7 SCSNA + 7 TAP = Q ^ 

Now we observe that in the case with Ui = (i = 1, . . . , N) the Onsager reaction term 
7 TAP Sj cancels out with the term 7 SCSNA Sj that emerges as a result of the evaluation of the 
correlation between the state of oscillators and the stored patterns. Then we assume Eq. fiffi) 
to hold generally so that ^ TOTAL = 0. As will be shown later the results obtained based on 
this assumption show good agreement with the results of numerical simulations. 



IV. BEHAVIORS OF THE NETWORK WITH A DISCRETE NATURAL 

FREQUENCY DISTRIBUTION 

For the sake of simplicity we focus on the behavior of the oscillator network with a 
discrete natural frequency distribution g{ui): 

where a represents the ratio of the number of oscillators with uji = to the total number of 
oscillators N. 



A. Appearance of a window for break down of the retrieval states 

To roughly sketch the effects of the natural frequency distribution with three frequency 
components, we investigate the behavior of the overlap with change of u)\ in the case of 
D = 0. In Fig. |3], we give c^-dependence of the overlap calculated from Eqs. (33)-(H0p and 



the result of numerical simulations with N = 4000 for the case with a = 0.02, a = 0.7, and 
D = 0. Good agreement between the theory and numerical simulations implies the validity 
of the present analysis. 

As is expected, an entrainment indeed occurs in the case of small u\ (c<Ji<0.5), where 
one has successful retrieval accompanying a large size of overlap. Even in the case of large 
Mi (LO^cji) successful retrieval can be realized with small size of overlap, since natural 
frequency of aN oscillators remains 0. 

To give an qualitative explanation for the occurrence of a window for break down of the 
retrieval states, we consider the case with a = for the moment. We can easily obtain the 
values of order parameters m, q, and u as a function of u>i (see appendix |TV| ) as shown in 
Fig. £|(a) . We see a phase transition to occur at u)\ ~ 0.7, and u is seen to increase as u)\ 
approaches u^: u — > 1 as U\ — >• while q = 1 for uj\ < io\ owing to the entrainment. 

Even when a ^ 0, such an enhancement of u around uj\ remains unchanged as can be 
seen in Fig. |](b). Noting Eq. (^) we can easily understand the noise variance ar/2 = 
aq/(l — u) 2 may be enhanced accordingly in the interval 0.5^co>i<1.0, where retrieval states 
disappears. 

In Fig. [| we draw the uj — a phase diagram to show the behavior of the storage capacity 
as a function of U\. The window observed in Fig. ^| turns to arise from the valley of a c (oui) 
curves. 
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B. Effect of white noise 



Behaviors of synchronization in the networks of coupled oscillators with white noise dif- 
fer from those of the deterministic networks with D = 0. In the absence of white noise 
(D = 0), one can in general divide the oscillators with Ui ^ into two groups of syn- 
chronized and desynchronized oscillators accordingly to the criterion of whether the phase 
velocity of an oscillation vanishes or not. Noting Eq. (36) and the form of effective transfer 
function with D = (Eq. ([K])) illustrated in Fig. |2|, we find that desynchronized oscilla- 
tors do not contribute to the value of overlap m, though they contribute to the value of 
other order parameters such as q and u. This is because we are concerned with the time 
averaged behavior of the local fields, where the time-averaged phase difference between the 
desynchronized oscillators with natural frequencies to and —u should be ir in the absence 
of white noise (i e. /(— u, h, h*) = —f(u, h, h*) for \h\ < uj). In the case with white noise 
(D ^ 0), however, the phase of each oscillator with ui ^ evolves with a certain non-zero 
time-averaged phase velocity, since the action of white noise prevents any oscillators from 
settling into fixed points. Hence it becomes impossible to distinguish between the synchro- 
nized and desynchronized oscillators any more. Nevertheless a finite value of the overlap is 
realized because of the existence of the equilibrium probability distribution on [0, 2ix) that 
is achieved after a long time with z ^ 0. 

In Fig. |^ we display the behavior of the overlap as a function of the noise intensity D 
obtained by the theoretical analysis and numerical simulations. As expected, the size of the 
overlap decreases as the noise intensity increases until the system undergoes a discontinuous 
transition at a critical noise intensity D c , above which a disordered state with m = is 
realized. Good agreement between the two results implies the validity of our treat ment 
based on the time averaged local fields together with the assumption that the TAP-like 
equation also holds in the case with a distribution of natural frequencies. 

In Fig. [7], we give U\ — a phase-diagram that represents the storage capacity plotted as 
a function of uo\ for various values of D. In most of the region of the storage capacity a c 
decreases as the noise intensity increases. We see that for D smaller than a certain critical 
D , a c (u>i) exhibits non-monotonic behavior with change in u>i, while for D > D a c (ui) is 
monotonically decreasing with U\. The occurrence of a window for the break down of the 
retrieval states with fixed a for D < D$ turns out to be attributed to the appearance of 
valley caused by the non-monotonic behavior of the a c (oj\) curve as in the case of D = 
(Fig-!). 



V. SUMMARY AND DISCUSSIONS 

We have investigated properties of an associative memory model of oscillator neural 
networks based on simple phase oscillators, where the influence of white noise together with 
a natural frequency distribution is considered in the case of an extensive number of stored 
patterns. In the presence of white noise every oscillator as well as its local field undergoes 
fluctuating motions even in the stationary state after a long time. To deal with such a 
situation we have taken an approach based on the concept of the TAP-like equation. To 
approximately derive the TAP-like equation for the system without an energy function we 
have taken the time average for the fluctuating local field of each oscillator neuron to make it 
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constant in time. On the basis of the time-averaged local field we have dealt with the single- 
body Fokker-Plank equation to obtain the time averaged outputs of the oscillators in the 
stationary state, from which we have evaluated the effective transfer function. The relation 
between the time-averaged output and the local fields involving such a transfer function can 
be viewed as an analogue of the naive TAP-like equation without considering the so called 
Onsager reaction term. We have supposed the proper form of TAP-like equation to be given 
by appropriately adjusting the Onsager reaction term such that setting Ui — (i — 1, . . . , N) 
naturally leads to the legitimate TAP equation, which we have obtained from the evaluation 
of the Gibbs free energy. Applying the SCSNA to this TAP-like equation, we have obtained 
the macroscopic order parameter equations, based on which the properties of associative 
memory of the network has been studied. 

Assuming a discrete symmetric natural frequency distribution with three frequency com- 
ponents for the sake of simplicity, we have presented the phase diagram showing the behavior 
of storage capacity as a function of the parameter uoi representing the width of the natural 
frequency distribution. In the case of D = the storage capacity a c has been found to 
exhibit non-monotonic behavior as U\ is varied and to attain a minimum at a certain u\. 
As a result of the occurrence of the valley in the uj\ — a c curve the break down of the re- 
trieval state with fixed a occurs for intermediate values of uji. When noise is present such 
a behavior has been found to be somewhat relaxed, and only for small values of the noise 
intensity D the phenomenon of the break down of the retrieval state can be observed. Our 
analytical result has shown excellent agreement with the results of numerical simulations. 

Our results show that associative memory based on temporal coding can be realized in 
the network of simple phase oscillators even in the presence of not only a distribution of 
natural frequencies but also external white noise. The result that temporal coding is robust 
against the existence of environmental noise is remarkable. Memory retrieval occurs in such 
a way that the oscillators undergo synchronized motions with the phase difference 0j — <pj 
between any two of the oscillators i and j settling into, for long times, the difference 0f — 9j 
of the memory pattern fi. In our model #f is chosen from uniformly distributed random 
numbers on [0,2tt). The resultant behavior, however, is qualitatively the same as that for 
our previous work (D = 0) on the special case where 6f = or 7r and hence is given 



by the well-known form = J2n€i£j with £f = ±1. A characteristic feature of memory 
retrieval accompanying synchronization is that, in contrast to networks with fixed point type 
attractors, each neuron exhibit oscillations in the local field or the membrane potential that 
are easily detected by other neurons in a certain network to determine whether memory 
retrieval is successful or not. Also worth noting is the appearance of two types of retrieval 
states with respect to the degree of synchronization: a high degree of synchronization that 
occurs for small u)\ with overlap m large and a low degree of synchronization that occurs for 
large u>\ with m small. 

The fundamental assumption we have used in the present study is the existence of 
the TAP-like equation ([16|) and ( |17D for our system together with the expression of 
7 TAP (Eq. (B3)). In the case of a;, = (i = 1,...,N) there occurs no problem because 
the genuine TAP equation exists as has been shown. In this case we have found that 7 TAP 
(Eq. (|I8D) is exactly canceled out by ^ SCSNA ; as in the case of the network of AGS, to yield 
7 TOTAL = together with the order parameter equations, that recover the ones by Cook [ 42] . 



who analyzed Q-state spin model including the case with Q — > oo for arbitrary temperatures 
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by means of the replica symmetric approximation. 

In the case with distributed natural frequencies, to obtain the form of the effective 
transfer function of the TAP-like equation, we replace time- dependent local fields hi in 
Eq. @ by their time-averaged ones, and assume that the Onsager reaction term of the 
form: 7 TAP Sj = — au/(l — u)si appears in effective local field as a result of fluctuation of 
local fields. The assumption of this form of the Onsager reaction term yields <-y TOTAL = 0, 
which leads to f(u, h, h*) = f(u, h, h*) in Eqs. 

Meanwhile the theoretical result we have obtained here is recovered by the TAP-like 
equation evaluated by simply replacing the local fields hi in Eq. (|5p by the time-dependent 
local field of the form: 

h . = h. + 7 SCSNA^ (46) 

where the renormalized local field: hi = £jm + (l/N) S M >i 12j£,i£,j* s j is constant in time. 
Note that, by use of this assumption, we can derive the order parameter equations without 
knowing the form of 7 TAP . In the present case both procedures give the same result, however, 
in some cases it seems that the later procedure gives the more accurate form of the transfer 
function. We will discuss details on this point somewhere else. 

Some special cases of the present model have also been investigated by several authors 
other than Cook. Arenas et al. |4l| have investigated the case with a = 0, where the natural 



frequency distribution is assumed to obey a Gaussian distribution. The result of this case 
can also be recovered by the present analysis. 

To our knowledge the case with distributed natural frequencies and a > was first 



studied by Park et al. |52] for different synaptic couplings by means of replica calculations 
based on the energy that is defined so as to satisfy 4>i — —dHp^/ dcpi + r\i. In our case 
the impark takes the form ifp ar k = — § Jij cos (4>i ~ 4>j ~ /%) ~~ Hi^i^i- However this 
energy does not make any sense because the equilibrium distribution exp (—H-p^/D) does 
not satisfy the periodical boundary condition P({(pi}) = P({<fii + 

Aonishi et al. |53| studied the case with D = and a Gaussian distribution for natural 



frequencies by considering that q = 1 holds in the set of SCSNA eqs. based on a different 
scheme from ours even in the presence of the group of the desynchronized oscillators. 

Yamana et al. also studied the deterministic oscillator network (D = 0) with a discrete 
distribution of natural frequencies that stores binary patterns by making an approximation 
that the motions of the group of desynchronized oscillators do not exert an influence on the 
behavior of the synchronized oscillators. Discarding the effect of desynchronized oscillators 
corresponds to considering the transfer function that takes the value zero inside the circle 
with radius u (see Fig. 0). For a wide class of natural frequency distributions this scheme 
seems to work to a good approximation in the case with D = 0, because the contribution 
of the desynchronized oscillators to such order parameters as m, q, and u is small. It is 
noted, however, that, in the case of D ^ 0, the phase of every oscillator with Ui ^ 
evolves with a certain non-zero time-averaged velocity and hence one cannot distinguish 
between synchronized and desynchronized oscillators. Accordingly for stochastic networks 
with D methods based on neglecting the effect of the desynchronized oscillators will not 
make sense and one has to deal with all of the oscillators equally as in the present analysis. 

Finally, we briefly discuss the relevance of our results to biologically related models of as- 



sociative memory. Biologically relevant models |^ , |32| , |54] -|57|j should be based on such spiking 



14 



neurons as the Ho dgkin- Huxley type [|54j] and integrate-and-fire type neurons |]31| , |32| , [56| . |57 . 



A simple integrate-and-fire neuron that is defined by 1-dimensional linear equation except 
for firing event can be described in terms of phase that is obtained by properly scaling the 
1-dimensional output variable. Synaptic couplings implimented into spiking neural networks 
are often assumed to incorporate the so called alpha function |56|] or its variant represented 
by the dynamics of a certain gating variable |5^j55[1 . So, major differences between the sim- 
ple phase oscillator model we have dealt with on the basis of the diffusive couplings among 
the oscillators and the spiking model will be the form of the synaptic couplings together 
with symmetry of an individual oscillator with respect to rotation of the phase variable. 
While the present model is assumed to take a sinusoidal phase interaction for simplicity, 
a spiking model with a synaptic interaction based on the alpha function takes the form of 
pulse like couplings pl| , p2| , |56| , p7| , which will lead to considering higher harmonics in the 



phase interaction. 

A spiking neural network model of associative memory we previously studied using 
FitzHugh-Nagumo neurons exhibits a nearly comparable size of the storage capacity to 
that of the standard analog network with the transfer function F(h) = (sgn(/i) + l)/2 that 
is larger than the storage capacity of the present model ||22|| . It will then be of interest to 
observe the outcome of introducing higher harmonics in the phase interaction of the simple 
phase oscillator model. We expect the storage capacity of the phase oscillator network to 
increase when the higher harmonics is taken into account. Such analysis is now under way. 

The problem of investigating properties of neurons synchronizing the envelope of a burst 
of spikes is also of interest, but is beyond the scope of the present paper, which aims studying 
the effects of such heterogeneities as a natural frequency distribution and external noise on 
the robustness of temporal coding in the oscillator network of associative memory. We 
consider that taking not only phase but also amplitude as variables for oscillatory neurons 
will provide a solvable model suitable for studying the case with such synchronization in 
networks of bursting neurons, which is also under way. 
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APPENDIX I: DERIVATION OF THE TAP EQUATION (|l6j) AND (|17|) IN THE 

CASE WITH uj = (I = 1, . . . , N) 

To obtain the TAP equation for the present model with the energy function (|2|) we follow 
the method of Plefka |]50| and Nakanishi [[[l[] used for the SK model and neural networks of 



Ising spins. 

The Hamiltonian (^) with a complex- valued external field Ri + Hi included reads 
H = aH — (Ri cos 4>i + h s ^ n 4>i) 

i 

= ~2^ , . Cij " ZiZj * ~ ? O^ 008 ^ + hsincpi), (1) 
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where a is introduced for the analysis below. Applying Legendre transformation to the free 
energy corresponding to the Hamiltonian H, one has 

G(a, { Sl }) = -/T 1 logTrexp (-fill) + E (R lXl + I iVi ), (2) 



where (3 = 1/D and Sj = Xi + iyi = (^coscf)^ +i(sm(f)ij = {zij . \- ■ ■/ denotes expectation 

with respect to the Hamiltonian H. 

We perform a Taylor expansion with respect to a 

G(a,{s l }) = Y,- r a n , (3) 

n=0 n - 

where = d n G/da n \ a=0 . Noting Ri + il i = dG/dx-i + idG/dyi, we rewrite Eq. (||) in the 
form 



G(a, {si}) 



-/T 1 log Tr exp { -a(3H + /? E E ( 

I n=0 i U - \ 

-/r 1 iogz + E -o— ^ + -5T- % 



cos ©j H — - — sin i 



dxi dyi 



-Mog^Trexp j/?E( 



-(3 1 log — Tr exp { (3 E ( cos & + sin <M f 

x exp { -a/?if + /3 E E — T cos 4>i + — sin 4>i \ 

+ h\n\\dx l X * + dy t Vl ) 
G (o) 

I \ a n ( dG^ dG^ 
-I3" 1 log (exp j -a/3# + /9 E E ( "o - cos & + ~q— sin ' 

I n=l i \ ^ 



= G<°> - /T' log (exp i-afiH + /3 £ ^A,J)„ (4) 

with 

^ = ^ E {(dGW/dxi + idG^/dyt) (zi - srf + (dGW/dXi + idG^/dyX{ Zi - Si )}, (5) 

where Z = Trexp (^r - cos0j + ^j^- sin^jj and (• • •) denotes expectation with 

respect to the Hamiltonian H with a = 0. 
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Noting Si = (^Zij = \Zi/ i from Eq. (Q), it follows 

1 

Then, from this equation and Eq. (S), one has 



H 



c ijSiSj . 



(6) 



G(a, {s,}) = G (0) + a(#) - /T 1 log (exp J a/3E + (3 ~r A 



n=2 



III 



where 



B 



5Z<+/(^« s i){ z j s j) ■ 



(7) 



(8) 



Evaluating G(a, {si}) by expanding this equation upto third order in a yields 
G(a, { Sl }) = + (H) Q a - ^( B \c? - ^(b 3 )/ + 0(a% 



Zi 5? 



where it is noted that \B/ Q = 

= ((^i - = 0, we have 

G(a, {s l }) 

1 



2 



(9) 



for every integer n > 1. Then, noting 



+ 



^ 2 {^(2, 0)^(0, 2)c* 2 + Ei(0, 2)Ej(2, 0), 

#3 



*2 

■Ji 



+2^(1,1)^(1,1)44} 



+ 



96 



£ {Ei(3, 0)E j (0, 3)c*f + Ei(0, 3)^(3, 0)c 



*3 
ji 



+3^(2, 1)^(1, 2)4 2 4 + 3^(1, 2)£,(2, 1)44 2 } 







^ £ {^(2, 0)^(1, 1)^(0,2)4^4 + ^(1, 1)^(2,0)^(0, 2) 



* * * 

C ik C ji C jk 



(ijk) 



+Ei(i, i)^(i, 1)^(1, i)444 + Ei(o, 2)^(2, 0)^(1, 1)444 



+^(2, o)^(o, 2)^(1, 1)444. + ^(i, i)^(i, i)E k (i, 1)44 



+^(1, 1)^(0, 2)E k (2, 0)444 + ^(0, 2)Ej(l, l)E k (2, 0)444} 
+Ofa 4 



(10) 



where Ei{n } m) = ^(^ — Sj) (4 — s*) y , and (ijk) denotes all combination to be taken 
so that either two of the indexes do not coincide (note that (ij) implies i 7^ j). Then 
substituting Eq. (|3]) into Eq. ( |10"1) yields, in the limit iV — > 00, 
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ni r u Mo) a v- * * aN(3(l - qf 2 aN/3 2 (l - q) 3 3 / 4 x 

G(a, {si}) = G u - - 2^ a — a +0[a ), (11) 

where q = jr J2i \ Si \ ■ Note that all the relevant terms higher than the term of first order in 
a under the limit N oo comes from the following in Eq. (|T0|) 

E E tl (l,l)E l2 (l,l)...E ln ^lK t2 c* 2t3 ...c* nll . (12) 

Since every higher order term than the first order one contains —a n f3 n ~ l (B n ^jJn\, one may 
expect that it yields terms of the form of Eq. fll2"D. Summarizing those terms we will have 

G(a, { Sl }) = G<°> - - E <M -aN^^—{i- qf a \ (13) 

2 ^ n=2 2 ™ 

Then, noting dq/dxi + idq/dyi = 2si/N, we obtain 

dG .dG 

dxi dyi 

dG<® .dG® _ TAP 
= _ + ,_-ag^- 7 «, (14) 

where 7 TAP = -aa{a/3(l - <?)/2}/{l - a/?(l - g)/2}. 
In the case of a = H becomes 

H = -E{( 5G(O V^)cos0 i + (dGW/dy^smfa). (15) 

i 

Thus we have 



hUjidGW/dxif + idGW/dyi 



IolPjidGW/dxtf + idGW/dyi 
OG^/dxi + idG^/dyi 



^{dG^/dxtf + {dGW/d yi f 
= /(o,<9G {0) /<9^ + idG (0) /dy h (dG^/dxi + z<9G {0) /%)*), (16) 

1 r 2 * 

where Ik{f) — — / exp (r cos (p) cos k(p dcp, and f(0,h,h*) is just the effective transfer 

V27T */o 

function we introduced in Eqs. (|^,(|T4"|), and (|15|). 

Considering the case with a = 1, from Eqs. ( [L4"D and (|T6|) we finally obtain the TAP 
equation: 
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Sl = f(0,hJ AP ,hJ AP *), (17) 

h TAP = ^ . + 7 TAP S . + R , + (lg) 

7 TAP = -a ^ l - q) ' 2 (19) 

7 a i-P(l-q)/2- ^ 



APPENDIX II: DERIVATION OF EQ. fl27|) 

To derive Eq. (P), we substitute Eq. © into Eq. (||) to obtain 

1 



» 1 + ^EE^X + 7^ ! 



M>1 j^i 



r jy2 Z_/ Vj Sj ^ 



M>1 jj^i k^j 



(>■■:■'■:') 



Sfc b k 



r jy2 Z-/ Z-/ M Vj Sj 



M>1 jV« fc^j 



(1) 



Utilizing the relations ^^j 1 * = 1, ^ Ya Cf = 0> an d so on > the fourth term of Eq. ([[]) becomes, 
in the limit N — > oo, 



iV2 



Z^ Z^ ^* ^ 



M>1 i^i fc^j 



A v v fVfv* „m J_ \ ^ ^/ 



uX 
iV 



(2) 



Following the almost same scheme the fifth term of the right hand side of Eq. (JJ) is shown 
to vanish in the limit iV — > oo. Substituting Eq. (|22|) into the left hand side of Eq. (JJ) we 
obtain Eq. ($2%. 



APPENDIX III: DERIVATION OF EQS. fl3|) AND fl39|) 



The Eq. (^) is straightforwardly derived from the definition of u by noting ( / || 



9{Re(h)} 6 9{lm(ft)} 
Eq. ( |32|) it is suffice to prove that w is real. 



and performing integration by parts. To show Eq. (|39j) from 
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To show u is real, note the rotationary symmetry structure of the form of trans- 
fer function (§) as is illustrated in Fig. ||. Because of this symmetry structure of 
fuo, h, h*) we also have f(oJ, re l6 , (re td ) J = e te f(uj, r, r) in the presence of non-zero complex 
7 TOTAL . One also immediately finds f(u,r,r) = f(—u,r,r)* and f(u,r,r) = f(—u,r,r) . 
Then it follows that f(u, re ie , (re w )*) = f(-co, (re 10 )*, re 10 )* and (x - iy)f(cu,h,h*) = 

|(x + iy)f(— to, h*, Jij | . On the other hand noting g(cJ) = g(—oj) and changing the vari- 
ables for integration, we have, from Eq. fl38|), 



aru 



Ux + iy)f(-u,h\h)\Y (1) 
Accordingly we have 

^°^ U = (({^ ~ ^)^(^'^'^*)} ^ = V^ u * ( 2 ) 



to conclude that u is real. 



APPENDIX IV: DERIVATION OF THE MACROSCOPIC ORDER PARAMETER 
EQUATIONS FOR THE CASE WITH D = AND a = 

In the case with D = and a = 0, substituting Eq. ( |T0| ) into Eq. fl2T|), we have 

!a < m < to\ 

a + (1 -a) v - 1 c^<m 

Using Eq. (fLOl) we also obtain, from Eqs. (|37| ) and (|39l), g and w as a function of m: 

(;=r+(1 _ a) (iii^Ei!) 2 o <m < Wl (2) 



) 2m 0<m<^x 



where we have noted w = //lle{|^}\\ = //^{^^-"(ff — r§£) w\ ^at * s obtained by 

representing the local field with the polar coordinate, i e. h — re ld . 

As u)\ approaches the point of phase transition from below, u increases as is shown 



in Fig. |](a) . At the phase transition point where m = a + (1 — a) ym 2 — uj\jm and 
^■|a + (1 — a)\jm 2 — uj\ /m j = 1, one has u = 1. 
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FIGURES 

FIG. 1. The typical time evolution of the local fields observed in numerical simula- 
tions with N = 8000, a = 0.02, D = 0. Natural frequencies are chosen so as to 
obey (a) a Gaussian distribution g(oj) = exp (— uj 2 /2a 2 ) /V^ira 2 with a = 0.3 and (b) 
g(u) = 0.15S(u + 1.4) + 0.76(w) + 0.15<5(w - 1.4). 

FIG. 2. Graphical representation of the effective transfer function f(u, h, 7j*^ in the case 
with D = and uj > 0. The output f(oj,h,h*^j is represented by a vector at the posi- 
tion h on the complex plane. In the region, where h < 1, oscillators get synchronized with 

f(uj,h,h*^j =1, while inside the circle oscillators get desynchronized with f(uj,h,h*^J < 1. In 
the case of u < the rotational direction of the flow pattern gets reversed owing to the property 

FIG. 3. wi-dependence of the overlap obtained from the present analysis is plotted together 
with the results of numerical simulations with N = 4000 for the case with a = 0.02, a = 0.7, D = 0. 
Since theoretical analysis is based on taking a time average of physical quantites of interest, the 



results of simulation are displayed in terms of time-averaged quantities \rnM\ = ii z i 



FIG. 4. (a) wi-dependences of the order parameters m, q, and u obtained from Eqs. (|l]),(Q), and 
(^|) are displayed in the case with a = 0.7 and a = 0. (b) same as (a) for the case with a = 0.02. 
The gap separating two types of retrieval state (the large u\ regimes and small u\ regimes) implies 
the disappearance of retrieval states. 

FIG. 5. <jj\ — a phase diagram representing the behavior of strage capacities for the various 
values of a in the case of D = 0. 

FIG. 6. D-dependence of the overlap m obtained from the present analysis is plotted together 
with the results of numerical simulations with N = 4000 in the case with a = 0.7, a = 0.01, uj\ = 0.3. 

FIG. 7. Same as in Fig. |5| for various values of D in the case of a = 0.7. 
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Figure 1: The typical time evolution of the local fields observed in numerical 
simulations with N = 8000, a = 0.02, D = 0. Natural frequencies are ch osen 
so as to obey (a) a Gaussian distribution g(co) = exp (— 0) 2 /2a 2 ) / V2no 2 with 
o = 0.3 and (b) g((o) = 0.158(©+ 1.4) + 0.78(a)) + 0.158(a)- 1.4). 
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Figure 2: Graphical representation of the effective transfer function 
/ ^00, h,hj in the case with D = and (0 > 0. The output / ^0), h,h* j is 
represented by a vector at the position h on the complex plane. In the region, 



where \h\ < 1, oscillators get synchronized with 



f((0,h,h 



1, while in- 
< 1. In the 



side the circle oscillators get desynchronized with 
case of G)< the rotational direction of the flow pattern gets reversed owing 
to the property / ^— CO, h,h*^j = j/ (a,7?,h\ j . 
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Figure 3: (0\ -dependence of the overlap obtained from the present analysis is 
plotted together with the results of numerical simulations with = 4000 for 
the case with a = 0.02, a = 0.7, D = 0. Since theoretical analysis is based on 
taking a time average of physical quantites of interest, the results of simulation 
are displayed in terms of time-averaged quantities. 
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Figure 4: (a) (Oi -dependences of the order parameters m,g, and u obtained 
from Eqs. (68), (69), and (70) are displayed in the case with a = 0.7 and a = 0. 
(b) same as (a) for the case with a = 0.02. The gap separating two types 
of retrieval state (the large (Oi regimes and small Cfli regimes) implies the 
disappearance of retrieval states. 
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Figure 5: G)i — a phase diagram representing the behavior of strage capacities 
for the various values of a in the case of D = 0. 
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Figure 6: D-dependence of the overlap m obtained from the present analysis 
is plotted together with the results of numerical simulations with Af = 4000 in 
the case with a = 0.7, a = 0.01, (Oi = 0.3. 
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